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We analyse small parameters in selected models of biological excitability, including 
Hodgkin-Huxley (1952) model of nerve axon, Noble (1962) model of heart Purkinje 
fibres, and Courtemanche et al. (1998) model of human atrial cells. Some of the 
small parameters are responsible for differences in the characteristic timescales of 
dynamic variables, as in the traditional singular perturbation approaches. Others 
appear in a way which makes the standard approaches inapplicable. Wc apply this 
analysis to study the behaviour of fronts of excitation waves in spatially-extended 
cardiac models. Suppressing the excitability of the tissue leads to a decrease in the 
propagation speed, but only to a certain limit; further suppression blocks active 
propagation and leads to a passive diffusive spread of voltage. Such a dissipation 
may happen if a front propagates into a tissue recovering after a previous wave, e.g. 
re-entry. A dissipated front does not recover even when the excitability restores. 
This has no analogy in FitzHugh-Nagumo model and its variants, where fronts can 
stop and then start again. In two spatial dimensions, dissipation accounts for break- 
ups and self-termination of re-entrant waves in excitable media with Courtemanche 
et al. (1998) kinetics. 
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1. Introduction 

The motivation of this study comes from a series of numerical simulations of spiral 
waves (Biktasheva, Biktashcv, Dawes, Holden, Saumarez & M.Savill 2003) in a 
model of human atrial tissue based on the excitation kinetics of Courtemanche, 
Ramirez & Nattel (1998) (CRN). The spiral waves in this model tend to break up 
into pieces and even spontaneously self-terminate (see fig. 1). 

No mathematical model of cardiac tissue is now considered ultimate or can 
claim absolute predictive power. The spontaneous self-termination may be rele- 
vant to human atrial tissue or may be an artefact of modelling. Understanding 
the mechanism of this behaviour in some simple terms would allow a more direct 
and certain verification. This is difhcult as the models are very complex and the 
events depicted in fig. 1 have many different aspects. Traditionally, such under- 
standing has been achieved in terms of simplified models, starting from axiomatic 
cellular automata description (Wiener & Rosenblueth 1946) through to simpli- 
fied PDE models (FitzHugh 1961, Nagumo, Arimoto & Yoshizawa 1962) which 
allow asymptotic study by means of singular perturbation techniques (Tyson & 
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Figure 1. Self-termination of a spiral wave in the CRN model. Red colour component: 
transmembrane voltage E, green colour component: gating variable Oi. Diffusion coefficient 
D = 0.03125mm'^/ms, preparation size 75 x 75mm. See also (Biktasheva et al. 2003). 



Keener 1988). This approach can describe some of the features observed in fig. 1, 
e.g. the "APD restitution slope 1" theory predicts when the stationary rotation of a 
spiral wave is unstable against alternans (Nolasco & Dahlen 1968, Karma, Lcvine & 
Zou 1994, Karma 1994). The relevance of the "slope-1" theory to particular models 
is debated (Cherry & Fenton 2004), but in any case it only predicts the instability 
of a spiral wave, not whether it will lead to complete self-termination of the spiral 
wave, its breakup, or just meandering of its tip. We need to understand how the 
propagation of a wave is blocked. This has unexpectedly turned out to be rather 
interesting. Some features of the propagation block in fig. 1 can never be explained 
within the standard FitzHugh-Nagumo approach. As this was the only well devel- 
oped asymptotic approach to excitable systems around, we had either to accept 
that this problem is too complicated to be understood in simplified terms, or to 
develop an alternative type of simplified model and corresponding asymptotics. 

We chose the latter. This paper summarizes our progress in this direction in the 
last few years (Biktashev 2002, Biktashev 2003, Biktasheva et al. 2003, Suckley & 
Biktashev 2003, Biktashev & Suckley 2004, Suckley 2004, Biktashev & Biktasheva 
2005). The results on the asymptotic structure of the CRN model are published for 
the first time. 



2. Tikhonov asymptotics 

The standard Tikhonov-Pontryagin singular perturbation theory (Tikhonov 1952, 
Pontryagin 1957) summarized in (Arnold, Afrajmovich, Il'yashenko & Shil'nikov 
1994) is usually formulated in terms of "fast-slow" systems in one of the two equiv- 
alent forms 

dX/dt = f(X,Y) dX/dT = ef(X,Y) 

edY/di = g(X,Y) dY/dT = g(X,Y) ^ ' ' 

where e > is small, t is the "slow time", T is the "fast time", t = eT, X is the 
vector of "slow variables" and Y is the vector of "fast variables". The theory is 
applicable if all relevant attractors of the fast subsystem are asymptotically stable 
fixed points. So the variables have to be explicitly classified as fast and slow, and 
the system should contain a small parameter which formally tends to zero although 
the original system is formulated for a particular value of that parameter, say 1. 

We consider Hodgkin & Huxley (1952) (henceforth referred to as HH) and Noble 
(1962) (henceforth N62) models. Both can be written in the same form, 

dE/At = -Y.i{E,h,m,n)/CM, 

Ay/dt = {y{E)~y)/Ty{E), y = m,h,n, (2.2) 
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Figure 2. Tikhonov asymptotics of HH. (a) Characteristic times n during an AP. Thin 
solid magenta line: the AP for a reference, (b) Phase portrait of the fast subsystem (2.4) 
at n = 0.37, h — 0.02. Solid red line: horizontal isocline. Dashed blue line: vertical iso- 
cline. Filled black circles: stable equilibria. Dash-dotted green line: stable separatrix of the 
saddle, the boundary of attraction basins. Black dotted lines: selected trajectories, (c) A 
three-dimensional view of the slow manifold (the surface) . Solid line: the fold line. Dotted 
line: the selected trajectory and its projections on coordinate walls, (d) AP in the original 
model, £ = 1, solid red line, and when the fast variables are made faster, e = 10"'', dashed 
blue line. See also (Suckley & Biktashev 2003). 



where 'Ei{E,h,m,n) = (gKn"^ + gKi{E)){E - Ek) + {gNam^h + gNai){E - ENa) + 
gi{E — El) is the total transmembrane current (/iA/cm^), t is time (ms), E is 
the transmembrane voltage (mV), Ek, k = Na,K,l are the reversal potentials of 
sodium, potassium and leakage currents respectively (mV), gk are the correspond- 
ing maximal specific conductances (mS/cm^), n, m, h are dimensionless "gating" 
variables, Cm is the specific membrane capacitance (/zF/cm^). y{E) are the gates' 
instantaneous equilibrium "quasi-stationary" values, and Ty{E) are the character- 
istic timescale coefficients of the gates dynamics (ms). Definition and comparison 
of parameters and functions used in (2.2) for HH and N62 models can be found in 
(Suckley & Biktashev 2003). 

To classify the dynamic variables by their speeds, we define empiric charac- 
teristic timescale coefficients, r^. For a system of differential equations dxi/dt = 

fi{xi, . . . ,xn), i ^ 1, . . . ,N, we define Ti{xi, . . . ,xn) = {dfi/dxi)'^ . The r's ob- 
tained for m, h and n in this way coincide with Tm,h.n in (2.2), and this definition 
can be extended to other variables, e.g. E in the case of (2.2). 



Hodgkin-Huxley. Fig. 2(a) shows how r's change during a typical action po- 
tential (AP) in the HH model. The speeds of E and m exchange places during the 
AP, as do the speeds of h and n, but at all times the pair {E, m) remains faster 
than the pair {h,n). This suggests introduction into system (2.2) of a parameter e 
which in the limit e ^ makes variables E and m much faster than n and h: 

edm/dt = {m{E) ~ m)/T„,{E), 

edE/dt = —Yji{E,h,m,n)/CM, 

dh/dt - {h{E)-h)/Th{E), 

dn/dt = {n{E) -n)/TniE). (2.3) 

The properties of this system in the limit e — > are shown in fig. 2(b-c). The 
fast transient; corresponding to the AP upstroke can be studied by changing the 
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independent time variable in (2.3) to T = t/ e and then considering the hmit e — > 
whieh gives the fast subsystem of two equations for m and E, 

dm/dT = {■m{E)-m)/Tra{E), 

dE/AT = -T.j{E,h,m,n)/CM, (2.4) 

in which the slow variables h and n are parameters as their variations during the 
onset are negligible. An example of a phase portrait of system (2.4) at selected 
values of h and n is shown in fig. 2(b). It is bistable, i.e it has two asymptotically 
stable equilibria, and a particular trajectory approaches one or the other depending 
on the initial conditions. The basins of attraction of the two equilibria are separated 
by the stable separatrices of a saddle point, which is the threshold between "all" 
and "none" responses. A fine adjustment of initial conditions at the threshold will 
cause the system to come to the saddle point. This is a mathematical representation 
of the excitation threshold in Tikhonov asymptotics. 

For different values of n and h, the location of the equilibria in the fast subsys- 
tem vary. All equilibria {E, m) at all values of n and h form a two-dimensional slow 
manifold in the four-dimensional phase space of (2.3) with coordinates {E, m, ft,, n). 
Projection of this two-dimensional manifold into the three dimensional space with 
coordinates {E,h,n) is depicted in fig. 2(c). It has a characteristic cubic folded 
shape, with two fragments of a positive slope (as it appears on the figure), sepa- 
rated by an "overhanging" fragment of a negative slope. The borders between the 
fragments are the fold lines, seen as nearly horizontal solid curves on the picture. 
The positive slope fragments consist of stable equilibria, and the negative slope 
fragment consists of unstable equilibria (saddle points) of the fast subsystem. 

The points of this manifold are steady-states if considered on the time scale 
T ~ 1 or equivalently i ~ e. On the time scale t ^ 1 these points are no longer 
steady states, but we observe a slow (compared to the initial transient) movement 
along this manifold, which explains its name. Asymptotically, the evolution on the 
scale t ^ 1 can be described by the limit e ^ in (2.3), which gives a system of 
two finite equations and two differential equations, 

= rn{E) — m, 

= Y,i{E,h,m,n), 

dh/dt = (h{E)-h)/Th{E), 

dn/dt = {n{E)-n)/T^{E), (2.5) 

The finite equations define the slow manifold and the differential equations define 
the movement along it. 

Fig. 2(c) shows a selected trajectory of system (2.3) corresponding to a typical 
AP solution. The only equilibrium of the full system (2.3), corresponding to the 
resting state, is at the lower, "diastolic" branch of the slow manifold. If the initial 
condition is in the basin of the upper branch, the trajectory starts with a fast initial 
transient, corresponding to the upstroke of the AP, then proceeds along the upper, 
"systolic" branch of the slow manifold, which corresponds to the plateau of the 
AP. When the trajectory reaches the fold line, a boundary of the systolic branch, 
the plateau stage is over. The fast subsystem is no longer bistable. The only stable 
equilibrium is now at the diastolic branch. So the trajectory jumps down. This is 
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Figure 3. Tikhonov asymptotics of N62. (a,b,d) Notations similar to fig. 2. (c) A 
two-dimensional view of the stable (black solid) and unstable (green dashed) branches 
of the slow manifold, and typical pacemaker potential trajectories (the limit cycles) : solid 
red (e2 ~ 1) and dashed blue (e2 ~ 10^'^), vertical dash-dotted line shows value n — 0.5 
for the fast phase portrait on (b). See also (Suckley & Biktashev 2003). 



repolarization from the AP and it happens at the time scale t ^ e. The trajectory 
then slowly proceeds along the diastolic branch towards the resting state. 

So, an inevitable feature of the asymptotics (2.3) is that the solution at smaller e 
has not only a faster upstroke, but also a faster repolarization, and the asymptotic 
limit of the AP shape is rectangular as opposed to the triangular shape in the 
exact model, see fig. 2(d). This is undesirable as it means that asymptotic formulae 
obtained in this way produce qualitatively inappropriate results. 

The practical importance of this excercise is limited as in HH the AP are not 
much longer than upstrokes. 



Noble 1962. This model is more relevant to cardiac AP. The speed analysis 
of N62 model, similar to the one we have done for HH model, reveals a different 
asymptotic structure but ultimately similar results. Fig. 3(a) demonstrates three 
rather than two different time scales. Variable m is the fastest of all, we call it 
"superfast". Of the remaining three, variables E and h are fast and variable n is 
slow. So we need two small parameters, ei to describe the difference between the 
fast and superfast time scales, and £2 to describe the difference between the slow 
and fast time scales. System (2.2) then takes the form 

e2eidm/dt = {rn{E) — m) / t„i{E) , 

e2dE/dt = —J^j{E,h,m,n)/CM, 

e2dh/dt - (h{E) - h)/ThiE), 

dn/dt = {n{E) -n)/Tn{E). (2.6) 

Consider first the limit ei — > 0. The superfast subsystem consists of one differential 
equation for m. It always has exactly one equilibrium which is always stable. So after 
a supcrshort transient, m{t) is always close to its quasi-stationary value rn(E{t)). 
Thus, with an error ~ ei we may approximate m by to and discard the first equation, 
i.e. adiabatically eliminate superfast variable to. 

With the remaining system of three differential equations for E, h and n, we 
consider the change of the time variable as before, t = e2T, and proceed to the limit 
£2 0. This produces the fast subsystem in the form 

dE/dT = -T.i{E,h,m{E),n)/CM, 

dh/dT = (h{E)~h)/ThiE). (2.7) 
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Fig. 3(b) shows a phase portrait of this system for a selected value of n when (2.7) 
is bistable. As there is one slow variable, all the equilibria of the fast subsystem 
form a one-dimensional manifold, i.e. a curve, in the three-dimensional phase space 
with coordinates {E,h,n). Its projection on the plane {h,E) is shown in fig. 3(c). 
Again the stable equilibria correspond to one slope (negative for the given choice 
of coordinates) and the opposite slope corresponds to the unstable equilibria. The 
branches are separated from each other by fold points. Again we have an upper, 
systolic branch, separated from the lower, diastolic branch, and the system has no 
alternative but to jump from one branch to the other in the time scale t ~ £2. As 
it happens, there are no true equilibria in the N62 model at standard parameter 
values, so these jumps happen periodically, producing pacemaker potential. This 
corresponds to the automaticity of cardiac Purkinje cells. Certain physiologically 
feasible changes of parameters may produce an asymptotically stable equilibrium 
at the diastolic branch, i.e turn an automatic Purkinje cell into an excitable cell. 

The systolic branch is separated from the diastolic branch, so in the asymptotic 
limit £2 ^ 0, if the upstrokes are fast, the repolarizations are similarly fast. This 
is in contradiction with the behaviour of the full model, which makes such an 
asymptotic analysis unsuitable. 

3. Non-Tikhonov asymptotics 

Noble 1962. To overcome the difficulties of the Tikhonov approach, we have de- 
veloped an alternative based on actual biophysics behind N62 and other cardiac 
excitation models, see fig. 4(a). The upstroke of an AP is much faster than the 
repolarization as the two processes are caused by different ionic currents. The up- 
stroke is very fast as the fast Na current causing it is very large. However all other 
currents, including the outward K currents that bring about repolarization are not 
as large, and there is no reason to tie the two classes of currents together in an 
asymptotic description. Mathematically, this means that in the right-hand side of 
the equation for E, only the term corresponding to the fast Na current is large and 
should have the coefficient £^^ in front of it. 

Next, the fast Na current is large only during AP upstroke, and remains small 
during other stages. The current is regulated by two gating variables, ni and h, 
and the quasi-stationary value of the specific conductivity of the current, defined 
as W{E) = rn^{E)h{E), is always much smaller than one. This happens because 
m,h € [0, 1] and rn^{E) is very small for E below a certain threshold voltage E^m 
and h{E) is very small for E above another threshold voltage Ey^ and Eh < Em- So 
the ranges of almost complete closure o{m^{E) and h{E) overlap. So whenever E 
changes so slow that m and h have enough time to approach their quasi-stationary 
values, the fast Na channels are mostly closed. The possibility for the opening of a 
large fraction of Na channels only exists during the fast upstroke, as m gates are 
much faster than h gates and have time to open before h close. 

Thus, the facts that m gate is much faster than h gate, rn'^{E) for E < Em, 
h{E) <C 1 for E > Eh and Eh < Em, are all related and reveal why the upstroke 
of the AP is much faster than all other stages. 

We adopt the hierarchy of times suggested by the formal speed analysis in the 
previous section, i.e. m is a superfast variable, E and h are equally fast variables 
and n is a slow variable. We keep the same notation £1 and £2 for the corresponding 
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Figure 4. Non-Tikhonov asymptotics of N62, excitable variant, (a) Functions of E illus- 
trating the small quantities taken into account, (b) Phase portrait of the fast subsystem, 
(c) Phase portrait of the slow subsystem, (d) Action potential solution of (3.2) for 62 = 1, 
solid red, and for £2 = 10~^, dashed blue. See also (Biktashev & Suckley 2004). 

small parameters. The small parameters should also ensure that m'{E) and h{E) 
are small m some ranges of E but not in others. We identify this smallness with £2 
rather than ei, as it is to compensate the large value of gtqa outside the upstroke, 
and the large value of gna is described by £2 . We denote the e2-dependent versions 
of m^{E) and h{E) as rfv'^E] 62) and h{E] £2). In this way we obtain 

eiesdm/dt = {m{E-e2)-m)/T,n{E), in{E-0)^M{E)e{E~E„,), 

CndE/dt = €^^gNa{ENa~E)rr?h + gK{EK-E)n^ 

+ gNa, (E) {ENa - E) + gK, [E) [Ek -E)+ gi {Ei - E) , 

e2dh/dt = {t{E;€2)~h)/Th{E), h{E;{)) ^ H{E)e{Eh - E), 

dn/dt = {niE)-n)/Tn{E), (3.1) 

where 6() is the Heaviside function. 

The limit ei ^ gives adiabatic elimination of the m gate, m « rn{E). 

The analysis of the limit e2 — > is complicated by a feature incidental to N62 
and not found in other cardiac excitation models. The small conductivity of the 
window current, W{E), is multiplied by a large factor g^a- In N62 the resulting 
window component of is comparable to other small currents and cannot be 
neglected outside the upstroke. The implications of this complication are analysed 
in (Biktashev & Suckley 2004). The result, in brief, is that the e2-dependent part 
of (3.1) can, in the limit £2 0, be replaced with a "modified N62" model: 



CndE/dt 
e2dh/dt 
dn/dt 



e^^gNaiENa ~ E) M^{E)e{E 
{H{E)e{Eh-E)-h)/Th{E), 
(n{E)-n)/T,,{E), 



Ern)h + gK [Ek - E) 



G{E) 
(3.2) 



where G{E) = gNa (ENa - E) W{E) + gNa^E) {ENa - E) + gK^E) {Ek - E) + 
gi {El - E), and as before M{E) w rn{E) for E > Em, H{E) « h{E) for E < Eh 
and W{E) ^m^{E)h{E). 

The limit 62 ^ of (3.2) in the fast time T = t/e2 gives the fast subsystem 



CMdE/dT 
dh/dT 



gNa {ENa ~ E) M^{E)e{E - Em)h 
{H{E)9{Eh-E)-h)/Th{E). 



(3.3) 



As intended, (3.3) takes into account only the fast sodium current and the gates 
controlling it, and everything else is a small perturbation on this timcscale. The 
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phase portrait of (3.3) is unusual, see fig. 4(b). The horizontal isocline (the red set) 
is not just a curve but contains a whole domain E < Em- The vertical isocline (the 
blue set) lies entirely within the red set, so the whole line h = H{E)0{Eh — E) con- 
sists of equilibria. An upstroke trajectory may end up in any of the equilibria above 
Erm so the height of an upstroke depends on initial conditions. For subthreshold 
initial condition, voltage remains unchanged in the fast time scale. Exactly what 
happens at the threshold E = iJ,„ depends on details of approximating function 
M{E), but in any case it does not involve any unstable equilibria. This is all dif- 
ferent from Tikhonov systems (see the paragraph after equation (2.4)) where the 
height of the upstroke is fixed, subthreshold potential decays in the fast time scale 
and the threshold consists of unstable equilibria and, if appropriate, their stable 
manifolds. In this sense, asymptotics of (3.3) give a new meaning to the notion of 
excitability, completely different from that in the Tikhonov systems. 

Let us consider the slow subsystem of (3.2). For any value of n we have a whole 
line of equilibria in the fast system h — H{E)9{Eh — E). The collection of such lines 
makes a two dimensional manifold in the three-dimensional space with coordinates 
{E,h,n). So the fast variable h can be adiabatically eliminated on the time scale 
t 1. Thus the slow subsystem, i.e. the limit £2 — > in (3.2), is 



The phase portrait of this system is shown in fig. 4(c). Further discussion of its 
properties can be found in (Biktashev & Sucklcy 2004). Notice that voltage E 
features in both the fast and slow subsystems, i.e. it is a fast or a slow variable 
depending on circumstances. This kind of behaviour is not allowed in Tikhonov 
asymptotic theory, so it is "a non-Tikhonov" variable. 

Courtemanche et al. 1998. CRN is a system of 21 ODE modelling electric 
excitation of human atrial cells, see (Courtemanche et al. 1998) for a description. 

Formal analysis of the time scales of dynamic variables by the same method as 
we used for HH and N62, reveals a complicated hierarchy of speeds, which changes 
during the course of the AP (see fig. 5(a)). From variables with smaller r's, we select 
those that remain close to their quasi-stationary values during an AP, and which 
can be replaced by those quasi-stationary values without significantly affecting the 
AP solution. We call them supefast variables. These include to, Ua and w. As before, 
we denote the associated small parameter ei. 

Next, we identify the fast variables with speeds comparable to the AP upstroke. 
This is also done by comparing the instantaneous values of the variables with the 
corresponding quasi-stationary values, and checking how their adiabatic elimination 
affects the AP, for the AP solution after the initial upstroke. In this way, we identify 
variables h, Oa and d as fast, with associated small parameter €2- 

Similar to N62, the transmembrane voltage is e^"'^-fast only during the AP up- 
stroke due to e^^-large values of /Na during that period, and is slow otherwise. This 
is due to nearly perfect switch behaviour of the gates to and h. The definition of 
/Na in this model is more complicated as there is also the j gate; however j is slow 
and does not change noticeably during the upstroke. 

These considerations lead to the following parameterization of the model: 
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Figure 5. Non-Tikhonov asymptotics of the Courtemanche et al. 1998 model, (a) Char- 
acteristic times during a typical AP potential solution, as indicated by the legend. Red 
solid line: te, blue dashed line: r's of superfast (m, Ua, w) and fast (h, Oa, d), green 
dash-dotted lines: r's of other, slow variables, (b) Phase portrait of the fast subsystem 
(3.6). (c) Action potential upstroke, ei = £2 = 1 (solid red), ei — 10"'^, £2 = 1 (dotted 
blue), ei = £2 = 10^^ (dashed green), in the fast time T = t/€2- (d) Same, for the whole 
AP in the slow time t. See also (Suckley 2004). 



CMdE/dt = 


- {e2^I-Ni,{E,m,hJ) + 


E'AE,...)), 


eie2dm/dt ~ 


{in{E;e2)-m)lT„,{E) 


m{E; 0) = M{E)0{E - Em 


e2dh/dt = 


{t{E; e2)-h)/Th{E). 


t{E;0)^H{E)e{Eh-E), 


eie2dua/dt = 


{u^{E)-Ua)/TuAE), 




eie2dw/dt = 


{w{E)-w)/t^{E), 




e2doa/dt = 


{o^{E)-Oa)/ToAE), 




e2dd/dt — 


(d{E) ~ d)/Td{E), 





(3.5) 

where I^^iE, m, h,j) = gNa'm^hj{E — Enu) is the fast Na current and 'E'j{E, . . . ) 
represents all other currents. Here we have shown only equations that contain ei or 
62- All other equations are the same as in the original model. 

As before, we adiabatically eliminate the superfast variables in the limit ei 0, 
and turn £2^0 in the fast time scale T = t/e. This gives the fast subsystem 

CndE/dT = gNaj{ENa~E)M^{E)9{E-Em)K 
dh/dT = {H{E)e{Eh-E)~h)/Th{E), 

dOa/dT = {oZiE)-Oa)/ToAE), 

dd/dT = (d{E)-d)/Td{E). (3.6) 

Notice that in (3.6) the equations for Oa and d depend on i?, but equations for 
E and h do not depend neither on Oa nor on d. Thus, the evolution equations 
for E and h form a closed subsystem. This system is identical to (3.3), up to the 
values of parameters, definitions of the functions of E and the presence of the slow 
variable j as the factor at the maximal conductance of the Na current, g^a- The 
phase portrait of (3.6), fig. 5(b), is similar to that of N62, fig. 4(b). So the peculiar 
features of the fast subsystem of N62 are not unique and are found in many cardiac 
models, including CRN. 

With a view of a practical application of approximation (3.6), it is interesting to 
test its quantitative accuracy. This is illustrated in fig. 5, panel (c) for the shape of 
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the upstroke in the fast time T ^ijt, and panel (d) for the shape of the AP in the 
slow time t. We see that the approximation of the AP is very good in both hmits 
ei — > and ti — > 0, except for the upstroke: e.g. the peak voltage is overestimated 
by about 13 mV. This is mostly due to the limit t\ — > 0, i.e. replacement of m with 
rn{E) . So the accuracy could be significantly and easily improved by retracting the 
limit £1 — > 0, which amounts to inclusion of the evolution equation for m instead 
of the finite equation m = rn{E). Most of the qualitative analysis remains valid. 
However, here for simplicity we stick to the less accurate but simpler case ei = 0. 
Notice that in fig. 5(b-d), we took M{E) = 1, H{E) = 1; the error introduced 
by that was small compared to other errors, particularly the error introduced by 
m ~ m{E). 

4. Application to spatially distributed systems 

Fronts. We now use the fast Na subsystem of the cardiac excitation (3.6) to 
consider a propagation of an excitation front through a cardiac fibre. In one spa- 
tial dimension, this requres replacement of ordinary time derivatives with partial 
derivatives and adding a diffusion term into the equation for E: 

dE/dt = J{E)e{E ~ E„,)h + Dd'^E/dx\ 

dh/dt = {e{Eh-E)-h)/ThiE), (4.1) 

where J{E) — QNajiENa - E)M^{E) and we have put H{E) — 1. We consider 
solutions in the form of propagating fronts. For definiteness, let us assume the fronts 
propagating leftwards, so E{x,t) = -E'(^), h{x,t) ~ h{^), = x + ct, h{—oo) = 
1, h{+oo) = 0, E{—oo) — E- and E{+oo) = £'+. In this formulation, we have 
three constants characterizing the solutions, the prefront voltage E^, the postfront 
voltage E+ and the speed c. It is not obvious which combinations of the three 
parameters admit how many front solutions. So we have considered a "caricature" 
of (4.1) by replacing functions J{E) and Th{E) in it with constants: 

dE/dt = J6{E - Er„)h + Dd'^E/dx^, 

dh/dt = {e{Eh-E)^h)/Th, (4.2) 

This system is piecewise linear and admits complete analytical investigation. 
Details can be found in (Biktashev 2002, Biktashev 2003); here we only briefiy 
outline the results. Fig. 6(a) illustrates a typical front solution. It exists if speed 
c and pre-front voltage E- satisfy a finite equation involving also the constants J 
and Th- The resulting dependence of the conduction velocity c on excitability J for 
a few selected values of E- is shown in fig. 6(b). These front solutions exist only for 
J at or above a certain minimum J,„i,i which depends on E^. For J > Jniin(£'-) 
there are two solutions with different speeds. Numerical simulations of PDF system 
(4.2) suggest that solutions with higher speeds are stable and solutions with lower 
speeds are unstable; this has been confirmed analytically by Hinch (2004). 

The replacement of functions J{E) and Th{E) with constants is a rather crude 
step. The purpose of the caricature is not to provide a good approximation, but to 
investigate qualitatively the structure of the solution set. To see if this structure is 
the same for the more realistic models, we have solved numerically the boundary- 
value problems for the front solutions in (4.1). There the role of the excitability 
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(a) (b) (c) (d) 

Figure 6. Fronts in the spatially distributed Na subsystem, (a) The structure of the 
front solution in the caricature model, (b) Speed of the fronts as function of excitability, 
at selected values of the pre- front voltage, in the caricature model (4.2). J* is minimal 
excitability at which propagation is possible at any E- , and c* is the corresponding prop- 
agation speed, (c) Speed of the fronts as function of excitability, at selected values of 
pre-front voltage, in the Na subsystem of the CRN model (4.1). On (b) and (c), solid 
red lines, above and raising, are the stable branches and dashed green lines, below and 
decreasing, are the unstable branches, (d) For comparison: speed of the fronts in a typical 
Tikhonov front (fast susbsystem of the FHN model). See also (Biktashev 2002, Biktashev 
& Biktasheva 2005). 

parameter is played by the variable j. The results of the calculations are shown 
in fig. 6(c). Not only the topology of the solution set is the same, but the overall 
behaviour of c(j, in (4.1) is quite similar to that of c(J, £■_) in (4.2), despite 
the crudeness of the caricature. 

PDE simulations show that approximation (4.1) overestimates the conduction 
velocity by almost 50% compared to the full model, and the error is again mainly 
due to the adiabatic elimination of the m gate. 

After eliminating the superfast variables m, Ua and w and the fast variables h, 
Oa and d, and retaining the non-Tikhonov variable _B, the slow subsystem of (3.5) 
has 16 equations. It describes the AP behind the front. 

The most important conclusion is that for any particular value of the prefront 
voltage E- there is a certain minimum excitability jmin = jminiE- ) and correspond- 
ing minimum propagation speed Cmin = Cmin (£<-), and for j < jmin, no steady front 
solutions are possible. This is completely different from the behaviour in FitzHugh- 
Nagumo (FHN) type systems, where local kinetics are Tikhonov and a front is a 
trigger wave in a bistable reaction-diffusion system. A typical dependence of the 
speed of such a trigger wave on a slow variable is shown in fig. 6(d): it can be slowed 
down to a halt or even reversed. The reversed trigger waves describe backs of prop- 
agating pulses in FHN systems. Thus, questions about the shape of the backs of 
APs and propagating pulses, and the spectrum of propagation speeds of a prop- 
agating excitation wave in a tissue come to be closely related. In both questions, 
our new non-Tikhonov approach provides different answers from the traditional 
Tikhonov/FHN approach. We have already seen that the new description is more 
in line with the detailed ionic models regarding the back of an AP. In the next 
section, we demonstrate the advantage of the new approach regarding the fronts. 

Dissipation of fronts. The fast subsystem of a typical spatially-dependent car- 
diac excitation model, discussed in the previous section, only provides part of the 
answer. This description should be completed with the description of the slow move- 
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Figure 7. Effects of temporary propagation blocks, (a) In the CRN model, (b) In the 
FHN model, (c) In the caricature Na front, (d) In the Na subsystem of the CRN model. 
See also (Biktashev & Biktasheva 2005). 



ment. The fronts are passing so quickly through every given point that the values 
of the slow variables at that point change little while it happens. Away from the 
fronts, the fast variables keep close to their quasi-stationary values. In our asymp- 
totics this means, in particular, that the fast Na channels are closed, and E is not 
a fast but a slow variable. Assuming absence of spatially sharp inhomogcncities of 
tissue properties, simple estimates show that outside fronts, the diffusive current 
is much smaller than ionic currents, so dynamics of cells there are essentially the 
same as dynamics of isolated cells outside AP upstrokes. 

Propagation of the next front depends on the transmembrane voltage E, which 
serves as parameter and the slow inactivation gate j of the fast Na current. 
This dependence gives an equation of motion for the front coordinate a;(i), 

iix/dt = c{j{x,t),E{x,t)), (4.3) 

where the instantaneous speed of the front c is determined by the values of E and 
j at the sites through which the front traverses (in case of E this is the value which 
would be there if not the front). This can only continue as long as the function 
c{j,E) remains defined, i.e while j(x,t) > ji„in(i?(a;, t)). If the front runs into a 
region where this is not satisfied, its propagation becomes unsustainable. 

What will happen then is illustrated in fig. 7(a), where the parameter g^a was 
varied in space and time. To make the effect more prominent, we did not use smooth 
variation, but put gNa = ^ in the left half of the interval for some time and then 
restored it to its normal value. The propagating front reached this region while it 
was in the inexcitable state. The result was that the sharp front ceased to exist, it 
"dissipated", and instead of an active front we observed a purely diffusive spread 
of the voltage. The excitability was restored a few milliseconds later, but the sharp 
front did not recover and diffusive spread of voltage continued, leading eventually 
to a complete decay of the wave. Note that the back of the propagating pulse was 
still very far when the impact that caused the front dissipation happened. 

This is completely different from the behaviour of a FHN system in similar 
circumstances, shown in fig. 6(b). There propagation was blocked for almost the 
whole duration of the AP. And yet when the block was removed, the propagation 
of the excitation wave resumed. Only if the block stays so long that the waveback 
reaches the block site and the "wavelength" reduces to zero, the wave would not re- 
sume. Such considerations have lead to a widespread, would-be obvious assumption 
that shrinking of the excitation wave to "zero length" is a necessary condition and 
therefore a "cause" of the block of propagation of excitation waves (Weiss, Chen, 
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Qu, Karagucuzian & Garfinkcl 2000). Comparison of panels (a) and (b) in fig. 6 
shows that this is far from true for ionic cardiac models, where such reduction to 
zero length happens, but only as a very distant consequence, rather than a cause, 
of the propagation block. The true reason for the block is the disappearance of the 
fast Na current at the front, observed phenomenologically as its dissipation. 

We expect that the condition j > jminiE) can also serve as a condition of 
propagation in the non-stationary situation on the slow space time/scale. Moreover, 
we conjecture that the dissipation of the front will happen where and when the 
front runs into a region with j < jminiE). This is illustrated by a simulation shown 
in fig. 6(c). It is a solution of the caricature system (4.2), where the excitability 
parameter J has been maintained slightly above the threshold Jinin(£'-) outside the 
block domain, and slightly below it within the block domain. As a result, the front 
propagation has been stopped and never resumed even after the block has been 
removed. A similar simulation for the quantitatively more accurate fast subsystem 
of the CRN model, (4.1), is shown on panel (d). Both agree with what happens in 
the full model on panel (a), and both confirm that the condition j < jmin is relevant 
for causing front dissipation. 



Break-ups and self-terminations of re-entrant waves In two spatial dimen- 
sions, the condition j < jmin(£') may happen to a piece of a wavefront rather than 
the whole of it. Then instead of a complete block we observe a local block and 
breakup of the excitation wave. This happened in the episode shown in fig. 8. 

The white dotted horizontal line on the top panels goes across the region where 
the propagating wave has been blocked and front has dissipated. The details of 
how it happened are analysed on the lower three rows, showing profiles of relevant 
variables along this dotted white line. The second row shows the profiles of E, which 
lose the sharp front gradient after t — 4100 ms. The third row shows the peaks of the 
spatial distribution of the product m^h; the sharpness of these peaks corresponds to 
the sharp localization of /Na at the front, and their decay accompanies the process 
of the front dissipation. The most instructive is evolution of the profile of the j 
variable shown on the bottom row. Consider the column t = 4100 ms. The gradient 
of j ahead of the front, i.e. to the left of the peak of m^h, is positive, and the front 
is moving leftwards, i.e. towards smaller values of j. That is, the front moves into 
a less excitable area, left there after the previous rotation of the spiral wave. To 
the right of the peak of /Na the gradient of j is negative which corresponds to the 
fact that j decreases during the plateau of the AP. Thus its maximal value at this 
particular time is observed at the front. This maximal value is, therefore, the value 
that should be considered in the condition of the dissipation, j < jmin{E). 

As soon as the front has dissipated {t « 4120ms), the profile of j starts to 
raise, so the maximum of the j profile observed at t = 4120 ms is the lowest one. 
From the fact that dissipation has started we conclude that this maximum is below 
the critical value jmin- Assuming that dissipation usually happens soon after the 
condition j < j^in is satisfied (simulations of (4.1) show that this happens within a 
few milliseconds), this smallest maximum value should be close to jmin, which gives 
an empirical method of determining jmin from numerical simulations of complete 
PDE models. For this particular episode the empirical value of jmin was found to be 
approximately 0.3. This is about 50% higher than the jmin predicted for the same 
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Figure 8. Analysis of a break-up of a re-entrant wave in a simulation similar to fig. 1. Top 
row: snapshots of the distribution of the transmembrane voltage, at the selected moments 
of time (designated above the panels). The other three rows: profiles of the key dynamic 
variables (designated on the left) along the dotted line shown on the top row panels, 
at the same moments of time. The scale of E is [— 100 mV, mV]. The scale of m^h is 
[0,0.15]. The scale of j is [0,0.5]. Cyan dash-dotted line on j panels represents jmin- See 
also (Biktashev & Biktasheva 2005). 



range of voltages by (4.1); we attribute this to the approximation m k, m which 
caused similar errors in the upstroke height and front propagation speeds. 



5. Conclusion 

Our new asymptotic approach for cardiac excitability equations has significant ad- 
vantages over the traditional approaches. The fast subsystem, represented by equa- 
tions (3.3) and (4.1), appears to be typical for cardiac models. This predicts that 
front propagation cannot happen at a speed slower than a certain minimum and 
at an excitability parameter lower than a certain minimum. When these conditions 
are violated the front dissipates and does not recover even after excitability is re- 
stored. We have obtained a condition for front dissipation in terms of an inequality 
involving prcfront values of j and E. This condition can be used for the analy- 
sis of break-up and self-termination of re-entrant waves in two and three spatial 
dimensions. 
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